Numerical approximations using Chebyshev polynomial expansions: El-gendi's 

method revisited 
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We present numerical solutions for differential equations by expanding the unknown function in 
terms of Chebyshev polynomials and solving a system of linear equations directly for the values of 
the function at the extrema (or zeros) of the Chebyshev polynomial of order N (El-gendi's method). 
The solutions are exact at these points, apart from round-off computer errors and the convergence 
of other numerical methods used in solving the linear system of equations. Applications to initial 
value problems in time-dependent quantum field theory, and second-order boundary value problems 
in fluid dynamics are presented. 
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I. INTRODUCTION 

A major problem in modern physics is to understand 
the time evolution of a quark-gluon plasma produced fol- 
lowing a relativistic heavy-ion collision. Although a mean 
field theoretical approach |l], ||, || can provide a reason- 
able picture of the phase diagram of quantum field the- 
ories, such studies do not include a rescattering mecha- 
nism, which would allow an out of equilibrium system to 
be driven back to equilibrium. As such, the past few years 
have witnessed a major effort concerning the search for 
approximation schemes jl], ^|. [7j which go beyond mean 
field theory. In the process, new numerical techniques 
were required in order to solve the ever more challenging 
systems of complex integro-differential equations. 

In this paper we revive and extend an old spectral 
method based on expanding the unknown function in 
terms of Chebyshev polynomials, which plays a cru- 
cial role in implementing our non-equilibrium field the- 
ory program. Finite-difference methods, though lead- 
ing to sparse matrices, are notoriously slowly conver- 
gent. Thus the need to use higher-order methods, like the 
nonuniform-grid Chebyshev polynomial methods, which 
belong to a class of spectral numerical methods. Then 
the resulting matrices are less sparse, but what is appar- 
ently lost in storage requirements, is regained in speed. 
We do in fact keep the storage needs moderate, as we can 
achieve very good accuracy with a moderate number of 
grid points. 

The Chebyshev polynomials of the first kind of de- 
gree n, T n (x) with n < N, satisfy discrete orthogonal- 
ity relationships on the grid of the extrema of T^(x). 
Based on this property, Clenshaw and Curtis || pro- 
posed almost forty years ago a quadrature scheme for 
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finding the integral of a non-singular function defined on 
a finite range, by expanding the integrand in a series of 
Chebyshev polynomials and integrating this series term 
by term. Bounds for the errors of the quadrature scheme 
have been discussed in |j| and reveal that by truncating 
the series at some order m < N the difference between 
the exact expansion and the truncated series can not be 
bigger than the sum of the neglected expansion coeffi- 
cients Jl0|]. This is a consequence of the fact that the 
Chebyshev polynomials are bounded between ±1, and 
if the expansion coefficients are rapidly decreasing, then 
the error is dominated by the m + 1 term of the series, 
and spreads out smoothly over the interval [—1, 1]. 

Based on the discrete orthogonality relationships of the 
Chebyshev polynomials, various methods for sol ving lin- 
ear and nonlinear ordinary differential equations Jfl[ and 
integral differential equations JT2| were devised at about 
the same time and were found to have considerable ad- 
vantage over finite-difference methods. Since then, these 
methods have become standard and are part of the larger 
family of spectral methods fl3|| . They rely on expand- 
ing out the unknown function in a large series of Cheby- 
shev polynomials, truncating this series, substituting the 
approximation in the actual equation, and determining 
equations for the coefficients. El-gendi |Q has argued 
however, that it is better to compute directly the values 
of the functions rather than the Chebyshev coefficients. 
The two approaches are formally equivalent in the sense 
that if we have the values of the function, the Chebyshev 
coefficients can be calculated. 

In this paper we use the discrete orthogonality relation- 
ships of the Chebyshev polynomials to discretize various 
continuous equations by reducing the study of the solu- 
tions to the Hilbert space of functions defined on the set 
of (N+l) extrema of T/v (x) , spanned by a discrete (N+l)- 
term Chebyshev polynomial basis. In our approach we 
follow closely the procedures outlined by El-gendi Q 
for the calculation of integrals, but extend his work to 
the calculation of derivatives. We also show that similar 
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procedures can be applied for a second grid given by the 
zeros of Tjy(x). 

In our presentation we shall leave out the technical de- 
tails of the physics problems, and shall refer the reader 
to the original literature instead. Also, even though our 
main interest regards the implementation of the Cheby- 
shev method for solving initial value problems, we present 
a perturbative approach for solving boundary value prob- 
lems, which may be of interest for fluid dynamics appli- 
cations. 

The paper is organized as follows: In Section II we re- 
view the basic properties of the Chebyshev polynomial 
and derive the general theoretical ingredients that allow 
us to discretize the various equations. The key element 
is the calculation of derivatives and integrals without ex- 
plicitly calculating the Chebyshev expansion coefficients. 
In Sections III and IV we apply the formalism to obtain 
numerical solutions of initial value and boundary value 
problems, respectively. We accompany the general pre- 
sentation with examples, and compare the solution ob- 
tained using the proposed Chebyshev method with the 
numerical solution obtained using the finite-difference 
method. Our conclusions are presented in Section |v|. 



II. METHOD OF CHEBYSHEV EXPANSION 

The Chebyshev polynomial of the first kind of degree n, 
T n (x), has n zeros in the interval [—1,1], which are lo- 
cated at the points 

'*■(*- h)\ 

x k = cos l— 2 -) , fc=l,2,...,n. (1) 



In the same interval the polynomial T n (x) has n + 1 ex- 
trema located at 



x^ = cos — , k — 0, 1, . . . , n . 



(2) 



The Chebyshev polynomials are orthogonal in the in- 
terval [—1,1] over a weight (1 — a; 2 ) -1 / 2 . In addition, 
the Chebyshev polynomials also satisfy discrete orthog- 
onality relationships. These correspond to the following 
choices of grids: 

• If Xk (k = 1,2,..., N) are the N zeros of Tn(x) 
given by (jj]), and iii,j<N, then 



N 



^ Ti(x k )Tj(x k ) = a l S lj , 

k=l 

where the constants cti are 

{ N , i = 0. 



(3) 



• If stk are defined by <M), then the discrete orthogo- 
nality relation is 



JY 



}/' Ti(xk)Tj(x k ) = 0i Sij , 

k=0 

where the constants ft are 
ft = 



(4) 



N 

N, i = 0,N. 



Here, the summation symbol with double primes denotes 
a sum with both the first and last terms halved. 

In general, we shall seek to approximate the values of 
the function / corresponding to a given discrete set of 
points like those given in Eos. dl], |J). Using the orthog- 
onality relationships, Eqs. (0, 0T), we have a procedure 
for finding the values of the unknown function (and any 
derivatives or anti-derivatives of it) at either the zeros 
or the local extrema of the Chebyshev polynomial of or- 
der N. 

A continuous and bounded variation function f(x) can 
be approximated in the interval [—1, 1] by either one of 
the two formulae 



or 



/(*) 



JV-l 
3=0 



N 



J2" b i T ^ ' 

3=0 



(5) 



(6) 



where the coefficients a,j and bj are defined as 



iV 



a i = mJ2 fMUxk) , j = 0, . . . , N - 1 , (7) 

(8) 



k=l 
N 



b i = nT," f@k)Tj{x k ), j = 0,. 



k=0 



and the summation symbol with one prime denotes a sum 
with the first term halved. The approximate formulae (|^) 
and (j^) are exact at x equal to Xk given by Eq. ([j]) , and 
at x equal to x k given by Eq. (||), respectively. 

Derivatives and integrals can be computed at the grid 
points by using the expansions (|| ||). The derivative 
f'(x) is approximated as 



N-l 



fix) « J2 n £ ' T * x *> T iW ' ( 9 ) 



and 



k=l 



N 



3=0 



N 



fix) « Y," /(^) n E" Tj{ik) T 'i {x) ■ (10) 



k=0 



3=0 
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Similarly, the integral f(t) dt can be approximated as 

/x n _ iV— 1 „a; 

/(t)d* « X! 7v E'^(^) / T i(*) d *> 

- 1 fc=l 3=0 

(11) 



/(*)d* « afE" 2 ^*) / ij(t)dt. 

1 fc=0 iV 3=0 

(12) 



3=0 



Thus, one can calculate integrals and derivatives based 
on the Chebyshev expansions (^) and @ , avoiding the 
direct computation of the Chebyshev coefficients (Q) or 
(0), respectively. In matrix format we have 



/(*) dt 



S [/] 



[/'(*)] « ^ [/] 
for the case of the grid (Q), and 



/(*) dt 
[/'(*)] 



£ [/] 



(13) 
(14) 

(15) 
(16) 



for the case of the grid (g), respectively. The elements 
of the column matrix [/] are given by either f(xk), k — 
1, . . . , N or f(xk), k — , . . . , TV. The right-hand side 
of Eqs. (p|, H) and (Q, |l|) give the values of the in- 
tegral J_ 1 f(t)dt and the derivative f'(x) at the cor- 
responding grid points, respectively. The actual values 
of the matrix elements Sij and Dij are readily available 
from Eqs. (§ 0), while the elements of the matrices S 1 
and D can be derived using Eqs. (|l0|, O). 



III. INITIAL VALUE PROBLEM 

El-gendi jl^] has extensively shown how Chebyshev ex- 
pansions can be used to solve linear integral equations, 
integro-differential equations, and ordinary differential 
equations on the grid (||) associated with the (N+l) ex- 
trema of the Chebyshev polynomial of degree N. Also, 
Delves and Mohamed have shown that El-gendi's 
method represents a modification of the Nystrom scheme 
when applied to solving Fredholm integral equations of 
the second kind. To summarize these results, we con- 
sider first the initial value problem corresponding to the 
second-order differential equation 

y"(x) + p(x)y'(x) + q(x)y(x) = r(x) , (17) 



with the initial conditions 



It is convenient to replace Eqs.(^) and ([l8|) by an in- 
tegral equation, obtained by integrating twice Eq. (|l7j ) 
and using the initial conditions (|18|) to choose the lower 
bounds of the integrals. Equations ( |l7| ) and (18) reduce 
to the integral equation in y(x) 



p{x )y(x )dx 



X fX 



y( x ) -yo-(x + i) \y' + p(-i)yo 

q(x") - p'(x")]y(x")dx"dx' 
r{x")dx"dx' , 



-l J-i 

x px' 



(19) 



which is very similar to a Volterra equation of the second 
kind. Using the techniques developed in the previous 
section to calculate integrals, the integral equation can 
be transformed into the linear system of equations 



,j -= 6 i:j + Sijp{xj) + [S ]ij q{xj)-p'{xj) 



N 



A if] = C, 
with matrices A and C given as 

A - *■■ ' d ■■->'- * 1 r£f2 
Q = g(xi) , i,j = 0, 1 

Here the function g{x) is defined by 

g(x) = Va + (x + 1) y' + p(-l)ya 
r(x") dx" dx' . 



(20) 



As a special case we can address the case of the integro- 
differential equation: 



y"(x) + p(x)y'(x) + q(x)y(x) 



K(x,t)y(t) dt, 
(21) 



with the initial conditions (|18|). We define the matrix L 
by 

Lij = SijK(xi,Xj), i,j = 0, 1, ... ,N. 

Then, the solution of the integro-differential Eq. (^l]) sub- 
ject to the initial values ( |l8| ) can be obtained by solving 
the system of N linear equations (po|), where the matrices 
A and C are now given by: 



Aij — Sij + Sijp{Xj) 

+ [S%j Q(xj) ~ p'{xj) 



[S 2 L] 



C, 



ya + {xi + 1) y' + p(- l)yo 



t/(-i) 



Vo 



y'(-i) 



y ■ 



(18) 



with ij = 0,1,... ,N. 

We will illustrate the above method using an example 
related to recent calculations of scattering effects in large 
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N expansion and Schwinger-Dyson equation applications 
to dynamics in quantum mechanics []l6| and quantum 
field theory |l7j , and compare with results obtained using 
traditional finite-difference methods. Without going into 
the details of those calculations, it suffices to say that the 
crucial step is solving an integral equation of the form 

G(t,f) = G (t,t')-2 f ne{Q(t,t")}G(t",t')dt" 



+2 / Q{t,t")Ke{G{t",t')}dt" , (22) 
Jo 

for G(t,t') at positive t and t'. Here, G(t,t'), G (t,t'), 
and Q(t,t') are complex functions, and the symbols IZe 
and Xm denote the real and imaginary part, respectively. 
In quantum physics applications, the unknown function 
G(t,t') plays the role of the two-point Green function in 
the Schwinger-Keldysh closed time path formalism p8| , 
and obeys the symmetry 



G(t,t') = - G(f,t) 



(23) 



where by G(t,t') we denote the complex conjugate of 
G(t,t'). Therefore the computation can be restricted to 
the domain t' <t. 

By separating the real and the imaginary parts of 
G(t,t'), Eq. ( p2[ ) is equivalent to the system of integral 
equations 



Ke{Go(t,t')} (24) 
-2 J Ke{Q(t,t")}1le{G(t",t')} dt" , 

(25) 



Tm{G (t,f)} 



-2 f TZe{Q(t,t")}Xm{G(t",t')} dt" 
Jo 



TZe{G(t,t')} 
Xm{G(t,t')} 



+2 / Xm{Q{t,t")}Ke{G{t" ,*')} dt" . 
Jo 

The first equation can be solved for the real part of 
G(t, t'), and the solution will be used to find Xm{G(t, t')} 
from the second equation. This also shows that whatever 
errors we make in computing lZe{G(t, t')} will worsen 
the accuracy of the Xm{G(t,t')} calculation, and thus, 
Xm{G(t,t')} is a priori more difficult to obtain. 

The finite-difference correspondent of Eq. ( p2| ) is given 

as 

f-i 

G(i,j) = G (i,j)-2J2 e k 1Ze{Q(i,k)}G(k,j) 



fe=i 



3-1 



+2^2e k Q(i,k)TZe{G(k,j)}, (26) 



fc=i 



where e k are the integration weights corresponding to the 
various integration methods on the grid. For instance, 
for the trapezoidal method, is equal to 1 everywhere 
except at the end points, where the weight is 1/2. Note 



that in deriving Eq. (|26j) , we have used the anti-symmetry 
of the real part of G(t,t') which gives lZe{G(t,t)} = 0. 

Correspondingly, when using the Chebyshev-expansion 
with the grid (|2j), the equivalent equation that needs to 
be solved is 



A'* 



G (i,j) - G(i,j) + 2j2S ik Ke{Q(i,k)}G(k,j) 

k=0 

N 

- 2^2§ jk Q(i,k)TZe{G(k,j)}. 

k=0 

In this case the unknown values of G(t,t') on the grid 
are obtained as the solution of a system of linear equa- 
tions. Moreover, the Chebyshev-expansion approach has 
the characteristics of a global method, one obtaining the 
values of the unknown function D(i,j) all at once, rather 
than stepping out the solution. 
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FIG. 1: TZe{G(t,t')} : Chebyshev/exact result (filled) versus 
the finite-difference result corresponding to the trapezoidal 
method (empty). 



FIG. 2: Xm{G(t,t')} : Chebyshev/exact result (filled) versus 
the finite-difference result corresponding to the trapezoidal 
method (empty). 

In Figs. and | we compare the exact result and the 
finite-difference result corresponding to the trapezoidal 
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method for a case when the problem has a closed-form 
solution. We choose 

Q(t,t') = -sin(i-f')+icos(t-i') , 
G (t,t') = (t - 1') cos(i - *') 

+i[cos(t - t') - (t + t') sin(t - t')} , 
G(t,t') = sin(t-t') +icos(i-t') . 

As we are interested only in the values of G(t,t') for 
t' < t, we depict the real and imaginary parts of G(t, t') 
as a function of the band index r = i{i — l)/2 + j, with 
j < i, used to store the lower half of the matrix. Given 
the domain < t < 6 and the same number of grid points 
(N=16), the result obtained using the Chebyshev ex- 
pansion approach can not be visually distinguished from 
the exact result, i.e. the absolute value of the error at 
each grid point is less than 10~ 5 . As expected we also 
see that the errors made on Xm{G(t, t')} by using the 
finite-difference method are a lot worse than the errors 
on lZe{G(t, t')}. As pointed out before, this is due to the 
fact that the equation for He{G(t, t')} is independent of 
any prior knowledge of Tm{G{t, t')} while we determine 
lm{G(t, t')} based on the approximation of 1Ze{G{t 1 1')}. 



IV. BOUNDARY VALUE PROBLEM 

In principle, the course of action taken in the previous 
section, namely converting a differential equation into an 
integral equation, also works in the context of a boundary 
value problem. Let us consider the second-order ordinary 
differential equation 

y"(x) + f[y'(x),y(x),x] = 0, x G [a,b] , (27) 

with the boundary conditions 

g[y(a),y'(a)]=c a , h[y(b),y'(b)} = c b . (28) 

No restriction on the actual form of the function 
f[y'(x),y(x),x] is implied, so both linear and nonlinear 
equations are included. 

We integrate Eq. (^) to obtain 



y '(x)-y'(a)+ / f[y'(x'),y(x'),x']dx' = 0. (29) 

J a 

A second integration gives 

y(x) - y(a) - (x - a)y'(a) 

f[y'(x"),y(x"),x"]dx"dx' = 0.(30) 



x rx 



The last equation is equivalent to Eq. (|19|). However, 
for an initial value problem, the values of y(a) and y'(a) 
are readily available. In order to introduce the bound- 
ary conditions for a boundary value problem, we must 
consider first a separate system of equations for y(a), 
y{b), y'{a) and y'(b), which is obtained by specializing 



Eqs. (29) and ( p0| ) for x = 6, together with the bound- 
ary conditions given in ( p8| ) . Then one can proceed with 
solving Eq. (^) using the techniques presented in the 
previous section. For instance, the Dirichlet problem 



y( a ) = y( b ) 

leads to the integral equation 



0. 

















a Ja 


la J 



x / f[y'{x"),y(x"),x"]dx"dx' = 0. 



(31) 



(32) 



Note, that one can also double the number of unknowns 
and solve Eqs. (p9h and (p5(]) simultaneously for y(x) and 

In this section however, we will discuss boundary value 
problems from the perspective of a perturbative ap- 
proach, where we start with an initial guess of the so- 
lution yo that satisfies the boundary conditions of the 
problem, and write y = yo + e, with e being a variation 
obeying null boundary conditions. We then solve for the 
perturbation e such that the boundary values remain un- 
changed. This approach allows us to treat linear and 
nonlinear problems on the same footing, and avoids the 
additional complications regarding boundary conditions. 

We assume that yo (x) is an approximation of the solu- 
tion y(x) satisfying the boundary conditions (|28|). Then 
we can write 

y( x ) = Vo[x) + e(x) , 

where the variation e(x) satisfies the boundary conditions 

g[e(a),e'(a)]=0, h[e(b), e'(b)] = . 

We now use the Taylor expansion of f[y'(x), y(x), x] 
about y{x) — yo(x) and keep only the linear terms in 
e(x) and e'(x) to obtain an equation for the variation 
e(x) 



e"(x) 



df[y'(x),y{x),x] 



dy'(x) 
df[y'(x),y(x),x] 



e'(x) 



dy(x) 



y{x)=yo(x) 
e{x) 



y(x)=y (x) 



= -y'o(x) - f[y' (x),y a (x),x] . 
Equation (^) is of the general form (|l7j ) 

e"(x) + p(x)e'(x) + q(x) e(x) — r{x) 

with 

df[y'(x),y(x),x] 



(33) 



p(x) 

q(x) 
r{x) 



dy'(x) 
9f[y'(x),y(x),x] 



dy(x) 



y{x)=y {x) 



y{x)=y {x) 



y'o(x) - f[y' {x),yo{x),x\ 
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Using the Chebyshev representation of the derivatives, 
Eqs. (H [l0|), and depending on the grid used, we solve 
a system of linear equations ( |20| ) for the perturbation 
function e(x). The elements of the matrices A and C are 
given as 

Aij = [D 2 ]ij + p(xi)Dij + q(xi)Sij, 
d = r( Xi ) , i,j = 1,2,... ,N, 

for the grid ([!]), and 



1.0 



.4 



[D%j + p(xi)Dij + q(xi)Sij, 
r(xi) , i,j = 1, ... ,N- 1 , 



for the grid (g). 

The iterative numerical procedure is straightforward: 
Starting out with an initial guess yo (x) we solve Eq. (|3^ ) 
for the variation e(x); then we calculate the new approx- 
imation of the solution 



0.6 





exact £ 




° 64 £ 




.... oo 0* 




- 16 i 


\ 









0.0 



FIG. 3: Numerical solutions of Eq. ( pq ) obtained for different 
values of N, using the Chebyshev expansion approach; we 



chose yo (x) = x , for — 1 < x < 1 



new 



,,old 



e{x) , 



(34) 



and repeat the procedure until the difference e(x) gets 
smaller than a certain e for all x at the grid points. 

It is interesting to notice that this approach can work 
even if the solution is not differentiable at every point 
of the interval where it is defined (Gibbs phenomenon), 
provided that the lateral derivatives are finite. As an 
example, let us consider the case of the equation 



xy'{x) - y{x) = 0, 



(35) 



which has the solution y(x) = \x\. In Fig. |^ we compare 
the numerical solutions for different values of N on the 
interval [—1, 1]. We see that for N = 64 the numerical 
solution can not be visually discerned from the exact so- 
lution. Eq. (|3^) is a good example of a situation when it 
is desirable to use an even, rather than an odd, number of 
grid points, in order to avoid any direct calculation at the 
place where the first derivative y'(x) is not continuous. 

We apply the perturbation approach outlined above 
to a couple of singular, nonlinear second-order bound- 
ary value problems arising in fluid dynamics. The first 
example 



y A (x) 



0, A > 0, 



(36) 



gives the Emden-Fowler equation when A is negative. In 
order to solve Eq. ([56]), we introduce the variation e(x) 
as a solution of the equation 



e"(x) 



A 



e(x) 



<Kx) 



The second example we consider is similar to a particular 
reduction of the Navier-Stokes equations EG] 



y"(x) - = 0. 

y*{x) 



(37) 



In this case, the variation e(x) is a solution of the equa- 
tion 



e"(x) 



0(g) 

Vo( x ) 



e'(x) 



- 2 wk v '° {x) e{x) 



In both cases we are seeking solutions y(x) on the interval 
[0, 1], corresponding to the boundary conditions 



2/(0) - 2/(1) 



. 



(38) 



Then, we choose yo{x) = sin(7ra;) as our initial approxi- 
mation of the solution. Given the boundary values (|3"g|), 
we see that the function f[y'(x),y(x),x] exhibits singu- 
larities at both ends of the interval [0, 1]. However, since 
the variation e(x) satisfies null boundary conditions, we 
avoid the calculation of any of the coefficients at the sin- 
gular points no matter which of the grids (|l|, ||) we choose. 
We consider the case when the above problems have the 
closed-form solution y(x) = x(l — x), with A = 1/2 in 
Eq. ( |36| ) . In Fig. |J we compare the exact result with the 
numerical solutions obtained using the Chebyshev expan- 
sion corresponding to the grid (|l|). 

The last example we consider arises in the study of 
ocean currents, specifically the mathematical explanation 
of the formation of currents like the Gulf Stream. Then, 
one has to solve a partial differential equation of the type 



dx 2 



dy 2 



a(x,y) 



d_ 

dx 



fay) = g(x,y) , (39) 



subject to null boundary conditions. To illustrate how 
the method works in two dimensions, we consider the case 
of a known solution u(x,y) — sin(7ra;) * sin(ny), defined 
on a square domain [0,1] x [0,1] with a{x,y) = 1, and 
compare the results obtained via a Chebyshev expansion 
versus the results obtained via a finite-difference tech- 
nique. We choose the function uq (x, y) = xy(l — x)(l— y) 
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FIG. 4: Chebyshev/exact solution of Eqs. (|(| and ( I3TI ) versus 
numerical solutions obtained using the Chebyshev expansion 
approach. 



as our initial guess. In Fig. ||we plot the exact result ver- 
sus the finite-difference result corresponding to the same 
number of points {n x — n y = N—8) for which the pro- 
posed Chebyshev expansion approach is not distinguish- 
able from the exact result. The number of iterations nec- 
essary to achieve the desired accuracy is very small (typ- 
ically one iteration is enough!), while the finite-difference 
results are obtained after 88 iterations. Of course, the 
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FIG. 5: Chebyshev/exact solution (filled) of Eq. ( J3S| ) versus 
the finite-difference result (empty) obtained for N = 8, as a 
function of the band index r = (i— 1)N + j. 

grid can be refined by using a larger number of mesh 
points. Then, the number of iterations increases linearly 
for the finite-difference method, while the number of it- 
erations necessary when using the Chebyshev expansion 
stays pretty much constant. In general, we do not expect 
that by using the Chebyshev expansion, we will always be 
able to obtain the desired result after only one iteration. 
However, the number of necessary iterations is compara- 
bly very small and does not depend dramatically on the 
number of grid points. This can be a considerable ad- 
vantage when we use a large number of grid points and 



want to keep the computation time to a minimum. 



V. CONCLUSIONS 

We have presented practical approaches to the numer- 
ical solutions of initial value and second-order bound- 
ary value problems defined on finite domains, based on 
a spectral method known as El-gendi's method. The 
method is quite general and has some special advantages 
for certain classes of problems. This method can be used 
also as an initial test to scout the character of the solu- 
tion. Failure of the Chebyshev expansion method pre- 
sented here should tell us that the solution we seek can 
not be represented as a polynomial of order N on the 
considered domain. 

The Chebyshev grids ([!]) and (J2j> provide equally ro- 
bust ways of discretizing a continuous problem, the 
grid (0) allowing one to avoid the calculation of functions 
at the ends of the interval, when the equations have sin- 
gularities at these points. The fact that the proposed 
grids are not uniform should not be considered by any 
means as a negative aspect of the method, since the grid 
can be refined as much as needed. The numerical solu- 
tion in between grid points can always be obtained by 
interpolation. The Chebyshev grids have the additional 
advantage of being optimal for the cubic spline interpo- 
lation scheme pl[ . 

The Chebyshev expansion provides a robust method of 
computing the integral and derivative of a non-singular 
function defined on a finite domain. For example, if both 
the solution of an initial value problem and its derivative 
are of interest, it is better to transform the differential 
equation into an integral equation and use the values of 
the function at the grid points to also compute the value 
of the derivative at these points. 
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FIG. 6: Standard deviation of the approximation for 
TZe{G(t, t')} as a function of the number of grid sites: Cheby- 
shev (filled) and finite-difference (empty) results. 

It is a well-known fact that spectral methods are more 
expensive that finite-differences for a given grid size, so 
in order to reach some specified accuracy there is always 
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FIG. 7: Standard deviation of the approximation for 
Tm{G(t, t')} as a function of the number of grid sites: Cheby- 
shev (filled) and finite-difference (empty) results. 
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FIG. 8: Standard deviation of the approximation for the solu- 
tion of Eq. (B9j) versus the number of mesh sites n x n y — N 2 : 
Chebyshev (filled) and finite-difference (empty) results. 



a tradeoff: finite differences need more points, but are 
cheaper per point. The goal is to reach a certain numeri- 
cal accuracy requirement as efficiently as possible. There- 
fore let us discuss here some computational cost consid- 
erations for sample calculations. We shall comment on 
two of the calculations presented in this paper: the cal- 
culation for the G(t,t') (see Eq. (p^)) and the solution of 
Eq. (H|). In Figs. ||, [7] and || we depict the convergence 
of the finite-difference and Chebyshev methods for ob- 
taining the approximate solutions, and we illustrate the 
elapsed computer time in Figs. [| and [lC]. The error of 
the Chebyshev expansion method decays exponentially 
as a function of N, while the error of the finite-difference 
method can be expressed as a power of N. For both 
methods, the running time depends exponentially of N. 
We conclude that indeed the execution time required by 
the spectral method increases faster with the number of 
grid points than the finite-difference method. However, 
in order to achieve a reasonable accuracy (e.g. a < 10 -6 ), 
the Chebyshev method requires only a small grid, and for 
this small number of grid points the computer time is ac- 
tually modest. All calculations where carried out on a 



(rather old) Pentium II 266 MHz PC. There were no ad- 
ditional numerical algorithms required for performing the 
finite-difference calculation, as these involve simple iter- 
ations of the initial guess. Due to the global character 
of the Chebyshev calculation, one needs to solve a sys- 
tem of linear equations. Since this is a sparse system of 
equations we have employed an iterative biconjugate gra- 
dient method jlO| for obtaining the numerical solution. 
For both problems, the sparsity of the relevant matrices 
is - 2/N. 
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FIG. 9: Execution time for obtaining an approximation for 
G(t, t') as a function of the number of grid sites: Chebyshev 
(filled) and finite-difference (empty) results. 
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FIG. 10: Execution time for obtaining an approximate solu- 



tion of Eq. (B9|) versus the number of mesh sites n x n v = 
Chebyshev (filled) and finite-difference (empty) results. 



N 



Most importantly, we have shown that the Chebyshev 
expansion is applicable to efficiently solving complex non- 
linear integral equations of the form encountered in a 
Schwinger-Dyson approach to determining the time evo- 
lution of the unequal time correlation functions of non- 
equilibrium quantum field theory. In this particular con- 
text, spectral methods have made possible for the first 
time to carry out complex dynamical calculations at next 
to leading order in quantum mechanics and field theory. 
Our results will form the basis for future studies of quan- 
tum phase transitions. 
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